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Preface 


The  Nuclear  Engineering  Test  Reactor  at  Vright- 
Patterson  Air  Force  Base  is  somewhat  unique  in  that  it 
was  designed  primarily  to  test  the  effects  of  high-in¬ 
tensity  radiation  on  materials.  Two  large  voids  (test- 
cells)  on  the  east  and  west  faces  of  the  core  provide  the 
means  of  conducting  such  experiments  on  very  large  items. 
The  calculation  of  the  flux  in  these  test-cells  is  hence 
of  the  utmost  importance.  1  decided  to  write  my  own  code 
for  this  calculation,  although  other  discrete  codes 
are  available,  since  the  principal  purpose  of  a  master's 
thesis  is  that  of  a  learning  device. 

When  this  project  was  first  begun,  I  had  planned  to 
calculate  the  flux  in  the  test-cells  of  the  NETR  in  three 
dimensions .  After  the  program  had  been  debugged  and  some 
preliminary  runs  made  in  two  dimensions,  it  was  found  that 
the  running  time  which  would  be  involved  was  too  expensive 
for  a  project  of  this  type.  Hence  the  remainder  of  my 
thesis  period  was  spent  in  trying  to  refine  the  program 
sufficiently  to  get  reasonable  results  in  xy-georaetry. 

I  wish  to  thank  Drs .  Carlson,  Lathrop,  and  Lee  of 
Los  Alamos  Scientific  Laboratory  for  their  kindness  in 
answering  my  Letters  concerning  difficulties  I  encounter¬ 
ed  in  the  theory;  Dr.  Bridgman  of  APIT  for  his  assistance 
and  criticisms;  and  Lt.  Max  Thompson  of  AFIT  who  executed 
the  PDQ  code  for  me.  Finally,  I  wish  to  express  my  spe¬ 
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cial  thanks  to  ay  wife  Marti  who  t3rp®d  this  thesis. 

Lee  H.  Livingston 
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Abstract 


This  report  describes  a  calculation  of  the  neutron 
flux  in  the  test-cells  of  the  Nuclear  Engineering  Test 
Reactor  at  Vright-Pattsrson  Air  Force  Base,  Ohio,  using 
the  discrete  approximation  to  the  transport  equation 
in  xy-geometry .  The  calculation  was  done  by  writing  a 
computer  program,  S4c40,  for  the  IBM  709^  using  three 
energy  groups  and  l600  mesh  points.  Fast  cross-sections 
were  generated  with  the  General  Atomics  code  GAM-I,  and 
thermial  cross-sections  were  calculated  by  hand  assuming 
a  Maxwell -Boltzmann  distribution.  The  results,  which  were 
compared  with  the  multigroup,  two-dimensional  diffusion 
theory  code  PDQ,  show  a  significantly  higher  thermal  flux 
over  the  entire  reactor  and  a  more  slowly  decreasing  flux 
for  all  groups  in  the  test-cell. 
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A  TRANSPORT  CALCULATION  OF  THE  FLUX 
IN  THE  NUCLEAR  ENGINEERING 
TEST  REACTOR  TEST-CELLS 

I .  Introduction 

This  paper  describes  a  calculation  of  the  neutron 
flux  in  the  test-cells  of  the  Nuclear  Engineering  Test 
Reactor  at  Wright-Pat terson  Air  Force  Base,  Ohio,  by 
means  of  the  S^  approximation  to  the  transport  equation. 
The  calculation  was  done  using  three  energy  groups  and 
l600  mesh  points  in  a  code  written  in  Fortran  for  the 
IBM  709^  computer. 

Statement  of  Problem 

The  NETR  is  a  light-water  moderated,  thermal  reactor 
with  MTR-type  fuel  elements.  On  the  east  and  west  sides 
of  the  core  are  two  large  voids  with  lead  and  steel  walls, 
surrounded  by  heavy  concrete  as  shown  in  Figure  1 .  These 
cells  are  approximately  ten  feet  deep,  seven  feet  high, 
and  seven  feet  wide. 

The  purpose  of  the  test-cells  is  to  provide  a  high- 
energy  neutron  flux  environment  for  irradiation  effects 
studies  of  various  materials.  The  cells  may  be  filled 
with  various  gases  and  are  equipped  with  small  gauge 
tracks  to  allow  massive  experimental  apparatus  to  be  in¬ 
serted.  A  cadmium  curtain  is  in  juxtaposition  with  the 
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core  to  reduce  the  thermal  flux  from  the  spectrum  enter¬ 
ing  the  cells.  In  the  region  of  the  cadoiius,  vhich  is  a 
strongly  absorbing  mediuta,  and  in  the  test-cells,  Vhere 
one  is  everywhere  less  than  one  mean  free  path  from  a 
boundary,  diffusion  theory  is  inapplicable;  hence,  trans¬ 
port  methods  should  be  *sed  to  calculate  the  neutron  flux- 

Hethod  of  Attack 

The  transport  equation  is,  in  principle,  capable  of 
predicting  the  neutron  flux  as  a  function  of  position  f , 
energy  E  t  and  direction  of  travel  .  Unfortunately, 
only  a  few  special  cases  can  be  solved  in  closed  form 
(Ref  6),  and  for  most  engineering  problems  one  must  turn 
to  some  means  of  approximating  the  solutions. 

Two  of  the  more  common  methods  are  the  (Ref  l4:33b) 
and  (Ref  9)  approximations.  These  methods  are  versa¬ 
tile  and  accurate  but  involve  a  great  deal  of  mathemat¬ 
ical  complexity.  The  discrete  S  approximation  gives  good 
results  (Ref  55539)  is  a  great  deal  simpler  in  form. 

It  can  be  programmed  easily  for  a  computer.  The  discrete 
difference  equations  are  developed  in  Chapter  II. 

Chapter  III  describes  the  calculation  of  the  constants 
in  the  problem  and  Chapter  IV  the  results  obtained.  The 
conclusions  drawn  from  this  work  and  recommendations  for 
future  work  are  given  in  Chapter  V. 
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II.  DerivatioK  of  the  Discrete  Approximation 

in  XY-Geoaeti~y 


The  derivation  of  the  discrete  S.,  approximation  to 
the  transport  equation  begins  with  the  steady-state  trans 
port  equation  for  monoenergetic  neutrons  (Refs  4,  16:219) 
of  energy  group  g»  which  is 


'  {>) 

where  N  =  H(x.  y,  fJi  .(j))  is  the  neutron  flux  and  S  =  S(x,y) 
is  the  source  term.  The  transport  equation  is  a  statement 
of  the  conservation  of  neutrons.  -One  equates  losses  due 
to  spatial  leakage  and  absorption  and  scatter¬ 
ing  collisions  (Zf  N^)  »  to  the  source  from  fission, 

external  sources,  and  scatter  into  the  group  and  direction 
under  consideration.  The  expansion  of  the  dot  product  of 
the  direction  vector  fi  with  the  gradient  operator  V 
rectangular  geometry  gives 


jsiN0COS0^iSIN0SlN0 


(2) 


If  one  definesCOS  0  as  jX  t  then  SIN 

One  now  operates  on  the  transport  equation  with  the  oper¬ 


ator 


Xi 


A 


^  I  j  pn  p 

. -'-TT-TT-  |ox  p  m  m  <  > 

xAyA^A^  I  J  J  I 


(3) 


to  create  a  difference  equation  in  the  flux  as  a  function 
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of  the  four  variables  x,  y,  yU.,  and  (j)  .  Consider  the 
first  term  of  the  resulting  equation 


Xj  Yj  Mn 

-L^  fix  IdY  fc  COS0  -^NlWM} 

/Ay^lAYUnUA  J  J  J  J 

^  ^  Yj-,  MnJi-- 


{'<) 


AxAyA^A0  j 


dY  /dM  W 


Yj-i 


/VR? 


where  the  assumption  has  been  made  that  the  integral  of 
the  partial  derivative  of  a  function  equals  the  change  in 
the  function  over  the  interval  of  integration.  The  inte¬ 
gration  over  the  other  three  variables  assumes  that  the 
integral  of  a  function  equals  the  function  evaluated  at 
an  average  value  of  the  variable  times  the  length  of  the 
interval.  Hence  one  may  write. 


Yj  Mn 

IdY  fdU  ( 

ti-i  '^n-i 


VhF<cos^4'^(X|,y,mA)  - 


(6) 
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where  the  quantity  ^  /  Indicates  an  average  value  of 

that  quantity  taken  over  the  range  of  integration.  The 
resultant  term  is 


<cos0g)(YFg> 

Ax  Ay 


A 


‘V 


(8) 


One  now  makes  the  approximation  that 


and 


<C0S(4)iC0S(| 


(9) 


n  Co) 

for  some  values  of  f)  and  (p^  ,  These  values  are  de- 

termined  later  in  the  chapter.  Thus  Eq  (8)  can  be  written 


as 


Y^[N!Xi.Yj^n4)-N|Xi-i'Yj.Mn-0f)] 
^  X 


(11) 


For  notationaL  simplification,  write  Eq  (ll)  as 

?|p|COS0j^[Nj"Nj_|]/A  where  the  absence  of  a  sub¬ 

script  implies  the  average  value  of  that  variable.  By 
a  similar  procedure,  the  second  terra  is  found  to  be 
^pSfN0£[iNj~!\lj_|]/A Y  .  on  the  third  term,  one 
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obtains 


AxAyA^A^ 


K:  Y; 

rl  /-J  /-'I  r 

dX  jdY  jdfijdp  [EtIx,y!  N(xX/t,0) 

^i-l  '*J-1  ^n-i^^-l 


=  ETlXi.YjlNiXi,Y;,jUn4)=^TNi 


(12) 


The  source  term,  S^,  is  composed  of  four  quantities; 
virgin  neutrons  from  fissions,  scattered  neutrons  from 
high-energy  groups,  "weak"  or  within-group  scatter,  and 
any  external  sources.  By  assuming  all  scatters  to  be  iso¬ 
tropic  this  source  term  is  dependent  only  on  position,  as 
was  stated  earlier.  This  assumption  is  reasonable  since 
fission  has  been  found  experimentally  to  be  very  nearly 
isotropic,  and  scattering  events  are  virtually  isotropic 
for  atoms  with  an  atomic  weight  of  10  or  greater.  There 
are  no  external  sources  in  this  problem.  Operating  on 
this  term  as  before,  one  has 


Xi  Y,  /«n 
dX  fdY  fd/U  /d(S|s(X,Yl 


Xi-|Yj-| 


=  S|XiYj)=S.  (,3) 


Eq  (1)  now  has  the  difference  form 


^nCOS0? 

Ax 


N?-N? 


.^nSIN0^ 


Ay 


‘>1 


(14) 
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Eq  (ill)  is  the  desired  difference  equation  which  re¬ 
lates  the  five  discrete  fluxes  as  shown  in  Figure  2.  One 
laay  obtain  difference  equations  relating  any  three  of  these 
discrete  fluxes  by  use  of  the  diamond  difference  technique. 


>  .  0 

< 

Nj 

A  Typical  Cell  in  the  Spatial  Grid 
Figure  2 


The  diamond  difference  assumption  is  that  2N,  =  N . +N ,  , 

11  1-1 

=  N .+N  (Ref  13568),  which  says  merely  that  the  average 
J  o 

of  the  flux  across  one  cell  in  any  ono  variable  is  the 
average  for  the  whole  cell. 

One  may  now  use  this  approximation  to  eliminate  two 
of  the  five  fluxes  in  Eq  (l4).  Eliminating  and 
one  obtains 
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which  may  be  written 


Alternately, 


one  may  eliminate 


and  N  .  , 

j-1 


as  follows 


Nj-Nj-Nj.|-+Nj]+^ZT[N|+Ni.|]  -S\  (17) 


which  may  be  written 


-T]nCOS0j^ 

.  Ax  Ay 


.  Ax  Ay 


27|nSIN0^' 


A 


Y 


Nj=Si 


(18) 


Eqs  (16)  and  ( 18)  are  the  two  working:  equations. 
These  working  equations  have  been  incorporated  into  a 
computer  program  called  S4c40  which  is  described  in 
Appendix  A. 

The  source  term  in  these  equations  is  calculated 
in  the  following  manner  (Ref  3s6).  From  an  inspection  of 
Eq  (1),  it  is  seen  that  N(x,  y,  jJi  t  (j)  )  is  the  flux  per 
steradian  with  the  units  (neutrons/cm  -sec-steradian). 
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The  isotropic  or  all-angle  riux  is 

Since  the  problem  is  assumed  to  be  symmetric  with  respect 
to  z,  only  one  hemisphere  need  be  considered  in  the  cal¬ 
culation.  The  areas  on  the  unit  sphere  about  each  dis¬ 
crete  direction  are  equal  (Ref  I),  and  hence  the  fluxes 
so  calculated  are  of  equal  importance  or  weight.  The  all¬ 
angle  flux  nCt)  is  thus 


lr)  =  27r 

fey.  "i  I'll 

“1 

].» 

_ 1 

(19) 


where  n  is  the  number  of  discrete  directions  taken.  The 
source  term  per  steradian  is 

n 

i=i 

For  an  Sg  calculation,  as  shown  in  Figure  3»  N=4,  and 
for  an  calculation,  as  shown  in  Figure  4,  Nt=12,  where 
N  is  the  number  of  discrete  directions  in  space  on  the 
hemisphere.  For  an  problem  in  xy-geometry,  Eq  (l4) 
may  be  written 


Ax 


Ni-N 


i-i 


+E7Nj  =  Sj  (2’) 


m=:i,2, — ,|i  i.ii.ia,— ,i2iii_7 
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These  discrete  directions  are  the  quadrature  points  for 
the  integration  Over  the  hemisphere.  (For  other  inte¬ 
gration  schemes  see  Refs  12,  13:^0,  and  155^36.)  If  one 
defines  a  flux 

n 

R'lf 

this  flux  can  be  put  directly  into  the  source  term.  Note 
that  »  differing  only  by  a  constant, 

and  hence  can  be  taken  as  a  relative  measure  of  the  all¬ 
angle  flux.  The  set  of  N'  (r")  for  the  mesh  of  discrete 
points  is  the  end  result  of  the  calculations. 

The  values  of  /Xpf^and  (j)^  are  chosen  so  that  equal 
areas  (solid  angles)  are  described  on  the  unit  sphere. 

It  remains  now  to  determine  the  values  of  !^XT\ 
be  used  in  Eqs  (l6)  and  (18).  This  is  done  by  demanding 
that  the  approximations  give  the  same  results  as 
diffusion  theory  in  regions  where  the  latter  is  appli¬ 
cable  (Ref  55538).  In  this  case  the  flux  is  isotropic, 
and,  neglecting  the  amplitude,  the  current  at  a  point  is 
given  by 


where  the  w^  are  the  point  weights  in  the  numerical  inte¬ 
gration  scheme  used.  The  /^pfi  can  be  written  (for  equal 
areas  on  the  unit  sphere,  as  shown  in  Figures  3  and  4)  as 

1 1 
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Mm=-i+^(m+i!,m=0,!,2,-^ 


(24) 


where 


Vi/Z 

~ — i 

m=| 


-  IIIIL+  i\X 

-  Z\2  ■  I  2 


(25) 


One  assumes  a  form  Mfyj"  2  ^  *  where  £  is  some 

constant  factor  to  be  determined.  The  weights  are 


(26) 


SO  that 


I-  ^  rrt^ 


(27) 


Q  '  4Q^I 

Substituting  the  expressions  for  w^  and  jUff)  into  Eq  (23) 
and  using  the  following  algebraic  identities 


7:2 


j 


E 

X=l 


x=-^j|i+il 


(28) 


X=i 


(29) 


and 


y  X^=i^j-(j+i)^(2j^+2j-i) 


(30) 
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cne  finds  that 

rv2  ^ 

^  (31) 

m=! 

The  solution  for  c  is 

^  V  n^-i?n 

Hence,  the  general  expression  for  /^pn  is 

''^n2+2n-2 

These  ’i^alues  are  tabulated  in  Table  IV  of  Ref  2:230. 

The  ^  are  not  weighted  but  taken  to  be  the  appropriate 
s:ean  values,  and  are  given  in  Table  V  of  Ref  2:231.  The 
values  for  this  problem  are  given  in  Table  I  below. 


Table  I 

Values  of  Angles  for  an  Calculation 
in  XY-Geomet^ 


m 

Mm 

Mm 

j 

1 

1 

0.8703883 

1 

TTA 

2 

37r/4 

2 

0.667 

0.3481553 

1 

TT/b 

2 

JTT/b 

3 

57T/8 

4 

iTlfB 
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111.  Calculational  Me thod 

As  in  all  reactor  calculations,  choosing  the  constants 
in  the  problem  and  developing  a  suitable  model  (the  grid 
Tor  the  numerical  calculation)  are  two  of  the  most  diffi¬ 
cult  decisions,  mainly  because  so  many  choices  are  possi¬ 
ble.  This  chapter  includes  the  model  used  for  the  NETR 
and  the  method  of  generating  cross-sections  for  ShChO. 

S4c40  Calculational  Grid 

A  sketch  of  the  NEJTH  is  shown  in  Figure  1 ,  Note  the 
presence  of  water  on  both  the  north  and  south  faces  of 
the  reactor,  of  the  test-cells  on  the  east  and  west  sides, 
and  of  the  cadmium  slabs  on  the  faces  of  the  core  nearest 
each  test-cell  to  reduce  the  theimial  fltuc.  Heavy  (barytes) 
concrete  surrounds  the  test-cells  for  biological  shield¬ 
ing.  As  Figure  5  illustrates,  these  features  have  been 
incorporated  in  the  grid.  It  should  be  noted  that  this 
problem  is  inherently  two-dimensional  in  character. 

The  final  grid  as  shown  in  Figure  6  was  based  on 
that  of  Garrett  and  Penaranda  (Ref  8:25)  with  some  modi¬ 
fication  as  noted.  Fewer  points  have  been  taken  in  the 
core  and  in  the  concrete,  and  the  aluminum  and  water 
regions  surrounding  the  core  have  been  homogenized.  Ap¬ 
proximately  the  same  number  of  points  has  been  used  in 
the  test-cell  (15  vs  17  in  the  x-direction) .  These 
Diodificatlons  were  made  to  keep  the  computational  time 
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to  a  minimum,  consistent  with  acceptable  accuracy. 

The  thermal  flux  in  the  interior  of  the  cadmuim  cur¬ 
tain  was  set  identically  equal  to  zero.  This  can  be  view¬ 
ed  as  an  additional  “boundary  condition"  imposed  on  the 
problem  and  was  done  primarily  to  suppress  oscillations 
in  that  region  and  to  aid  convergence.  Carlson  predicts 
(Ref  2:222)  that  when  the  product  of  the  cross-section 
and  the  spatial  increment  is  large,  oscillations  in  the 
calculated  fluxes  may  occur.  This  was  found  to  happen 
both  in  the  cadmiiun,  and  in  the  more  removed  parts  of  the 
test-cell  when  an  increment  of  37  centimeters  was  taken. 
When  this  increment  was  reduced  to  22.5  centimeters  the 
values  were  more  stable,  although  some  oscillation  re¬ 
mained  . 

Cross-sections 

The  fast  and  epithermal  cross-sections  used  in  this 
calculation  were  calculated  with  the  GAM-I  code  (Refs  10, 

1 1 )  and  the  thermal  cross-sections  were  computed  by  hand 
assximing  a  Maxwell-Boltzmann  distribution  and  a  tempera¬ 
ture  of  58.5°  C.  (See  Appendix  B.) 

The  energy  limits  were  taken  to  be  the  same  as  those 
used  by  Garrett  and  PeTEaranda  (Ref  O)  for  purposes  of 
comparison  with  PDQ  (Ref  1 :78l-783»  787)»  Atom  densities 
of  the  various  materials  were  taken  from  manufacturer’s 
datai  the  density  of  air  was  computed  at  a  temperature 
of  20°  C. 
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Input  Flux 

The  initial  estimate  of  the  flux  was  calculated 
using  the  S2  approximation  and  a  20x20  raesh.  After  the 
S4c40  code  was  first  run,  the  last  computed  values  of 
flux  were  used  as  input  for  the  succeeding  run.  This 
allowed  the  problem  to  be  stopped  and  restarted  at  will. 
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IV,  Results 

The  sophistication  of  the  equations  in  the  last  chap¬ 
ter  and  the  complexity  of  their  couplin^f  required  that  a 
great  deal  of  effort  and  time  (about  2^  months)  be  spent 
in  debugging  the  computer  code.  The  correspondence  in 
Appendix  C  refers  to  five  of  the  many  problems  which  arose 
in  this  work.  This  correspondence  has  been  included  in 
the  hope  that  it  may  prove  useful  to  others  in  future 
work.  A  test  problem,  which  is  described  below,  was  em¬ 
ployed  during  this  debugging  to  gauge  the  output  of  the 
program. 

Test  Problem 

The  purpose  of  the  test  problem  was  to  determine 
whether  the  errors  in  Fortran  programming  had  been  elim¬ 
inated  and  to  judge  whether  the  output  of  the  program 
was  reasonable.  The  test  problem  was  an  approximation, 
using  first  a  water-U^^^  core  surrounded  completely  by 
water,  and  then  with  an  air  gap  in  the  water  reflector 
parallel  to  one  face  of  the  core  as  shown  in  Figure  7* 

The  atom  density  of  fuel  and  water  in  the  core  was  cal¬ 
culated  using  modified  one-group  theory.  Various  grid 
spacings  were  employed,  the  final  one  being  ten  points 
in  the  y-direction  10  cm.  apart,  and  sixteen  points  in 
the  x-direction  6.25  cm.  apart.  The  Initial  flux  esti¬ 
mate  was  made  intentionally  poor  to  determine  if  the 
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program  were  sensitive  to  these  input  values.  It  was 
found  that  the  initial  flux  strongly  influenced  the  cal¬ 
culated  values,  the  flux  shape  changing  very  slowly.  One- 
fourth  the  reactor  was  considered,  using  symmetry  along 
the  centerlines  in  both  the  x-  and  y-directions . 

The  calculated  flux  along  the  x-centerline  is  shown 
in  Figure  8.  These  values  were  computed  after  176  iter¬ 
ations.  The  thermal  flux  is  somewhat  high  in  the  core, 
being  approximately  equal  to  the  fast  flux.  The  values 
in  the  air  gap  look  very  good,  the  higher-energy  fluxes 
dropping  almost  exponentially  and  the  thermal  flux  drop¬ 
ping  initially  and  then  rising  because  of  backscatter  and 
thermalization .  A  two-group,  one  dimensional  diffusion 
theory  hand  calculation  gave  nearly  the  same  results. 

NETR  Results 

The  final  flux  calculated  along  the  centerline  of 
the  test-cell  is  shown  in  Figure  9;  for  comparison,  the 
flux  calculated  by  PDQ  is  shown  in  Figure  10.  (The  small 
insert  in  Figures  9-13  shows  the  location  in  the  test¬ 
cell  of  the  abscissa  of  the  figure.)  The  fast  and  epi¬ 
thermal  fluxes  are  in  fairly  good  agreement,  although 
the  epithermal  flux  is  noticeably  lower  in  magnitude. 

The  thermal  flux  predicted  by  SkCkO  is  two  orders  of 
magnitude  greater  than  that  predicted  by  PDQ;  this  is 
true  not  only  in  the  test-cell  but  also  for  the  entire 
reactor.  The  thermal  flux  was  everywhere  decreased  by  a 
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Relativ-'  Fluxes 
along  Test-Cell  Centerline 
(PDQ) 

Figure  10 
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Tactor  or  100  and  ShCkO  executed  again  to  see  If  tills  vere 
the  result  of  a  poor  initial  estimate.  After  five  iter¬ 
ations,  the  program  calculated  essentially  the  original 
values.  Due  to  time  limitations,  the  solution  of  this 
problem  must  be  left  to  other  investigators. 

ShekO  predicts  a  flux  decrease  of  one  order  of  magni¬ 
tude  in  100  centimeters  in  the  test-cell,  vhereas  PDQ 
predicts  a  drop  of  three  orders  of  magnitude  in  the  same 
distance.  This  indicates  that  a  higher  flux  i«  available 
for  experiments  than  had  been  thought  previously. 

The  flux  across  the  test-cell,  for  tvo  distances  in 
the  x-direction,  is  shovn  in  Figures  11  and  12.  The  re¬ 
lative  magnitudes  of  the  different  energy  fluxes  remains 
essentially  constant.  Near  the  neutron  vindow  the  flux 
drops  off  slowly  as  one  moves  away  from  the  centerline. 

At  a  distance  of  6j,6  centimeters  the  reverse  is  true, 
a  minimum  occuring  near  the  centerline  nearly  two  orders 
of  magnitude  lower  than  near  the  sides  of  the  test-cell. 
This  point  is  discussed  further  in  Chapter  V,  The  dip 
in  the  flux  shown  in  Figure  12  at  about  31  centimeters 
is  probably  due  to  shadow-shielding  by  the  test-cell 
wall.  (See  Figure  5*)  The  semi-logarithmic  plot  in 
Figures  11  and  12  is  misleading.  The  derivative  of  the 
flux  at  the  centerline  in  both  cases  is  actually  zero. 
Figure  11  has  been  plotted  on  a  linear  scale  in  Figure 
13  to  demonstrate  this  point. 
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V.  Conclusions  and  Recomiasndationa 

The  discussion  in  Chapter  XV  suggests  that  the  trans¬ 
port  equation,  although  basically  the  proper  approach  to 
the  problem,  predicts  fluxes  which  are,  in  this  case,  dubi¬ 
ous  at  best.  If  the  fliax  predicted  by  the  program  is  cor¬ 
rect,  one  possible  explanation  for  the  high  thermal  fliax 
in  the  test-cell  is  backscatter.  The  high-energy  neutrons 
emerging  from  the  core  may  become  thermalized  in  the  water 
and  concrete  and  subsequently  diffuse  back  into  the  test- 
cell,  thus  defeating  the  purpose  of  the  cadmium  curtain. 
This  hypothesis  could  be  tested  by  executing  S4c40  with 
a  cadmium  lining  around  the  test-cell,  (Appendix  A) 

The  S  approximation  is  easily  programmed  but  is 

Xl» 

rather  slow  to  converge,  Carlson  predicts  that  about  200 
Iterations  are  required  for  convergence  for  a  thermal 
reactor  (Ref  5^538)?  since  S4C40  requires  one  minute  six 
seconds  per  iteration,  this  would  mean  about  three  hours 
of  computing  time  to  achieve  convergence. 

There  are  a  number  of  refinements  which  may  be  in¬ 
cluded  in  the  program  to  increase  the  accuracy  of  the  re- 
sidts.  Four  will  be  mentioned  briefly.  The  most  signi¬ 
ficant  improvement  would  be  obtaining  more  accurate 
cross-sections.  In  the  core  region,  a  separate  trans¬ 
port  calculation  could  be  performed  to  predict  fission 
cross-sections  which  take  into  account  the  thermal  flux 
depression  in  the  fuel.  Total  and  scattering  cross- 
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sections  in  the  core  could  be  similarly  calculated,  per¬ 
mitting  the  homogenized  core  to  better  reflect  reality < 

In  addition  it  should  be  advantageous  to  divide  the  energy 
spectrum  into  more  groups.  Table  II  gives  a  suggested  set 
of  energy  groups  over  the  range  of  10  Mev  to  0.4l4  ev. 

These  groups  were  selected  to  reflect  the  resonance  regions 
of  the  predominant  elements  present  in  the  reactor.  It 
may  also  be  necessary  to  consider  group  structure  within 
the  thermal  region.  Such  data  could  be  generated  by  a 
program  such  as  SLOP-1  (Ref  1:7^5)* 

A  second  refinement  would  be  to  consider  anisotropic 
scatter,  especially  in  the  test-cell  where  the  relatively 
light  elements  oxygen  and  nitrogen  are  present.  Both 
Carlson  (Ref  3;19)  and  Lee  (Ref  13^110)  discuss  techniques 
of  handling  this. 

When  S4C40  calculates  a  negative  angular  flux,  the 
flux  is  set  to  zero  and  the  computation  continued.  A 
better  technique  would  be  to  adjust  the  angular  flux  on 
the  opposite  edge  of  the  cell  and  the  average  flux,  as 
well  as  setting  the  negative  quantity  to  zero  (Ref  4:37)* 

Still  a  fourth  refinement  would  be  to  adjust  the 
model  in  any  of  several  ways.  Since  the  NETR  is  not  ex¬ 
actly  symmetric,  as  shown  in  Figure  1,  a  model  using  half 
the  reactor  would  be  more  realistic  (Ref  8:26,27).  Al¬ 
ternately,  more  mesh  points  in  the  -^-core  model,  which  is 
shown  in  Figure  6,  or  a  three-dimensional  model  using  five 
or  so  points  in  the  z-direction  would  improve  accuracy. 
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Table  II 

Cross-section  Group  Structure 


Group 

Lower  Energy 
(ev) 

Au 

Lower  GAM 
Group 

1 

10^-7.79x10^ 

0.25 

1 

2 

3.68x10^ 

0.75 

4 

3 

4.98x10^ 

2.00 

12 

k 

3.02x10^ 

0.50 

14 

5 

3.18x10^ 

2.25 

23 

6 

9.61x10^ 

3.50 

37 

7 

4.54x10^ 

0.75 

40 

8 

1.37x10 

3.50 

54 

9 

5.04 

1 .00 

58 

10 

0.4l 

2.50 

68 
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Any  of  the  refinements  listed  above  should  improve 
the  accuracy  of  the  code  and,  with  sufficient  iterations, 
any  desired  degree  of  precision  could  be  achieved. 
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Appendix  A 
The  S4CU0  Code 


General  Comments 

The  S4c40  code  was  written  Tor  a  particular  engineer¬ 
ing  application  and  was  not  designed  for  use  with  any  other 
problems.  It  can  be  modified  easily,  however,  to  incor¬ 
porate  the  suggestions  mentioned  in  Chapter  V.  The  pro¬ 
gram  is  written  to  accomodate  three  energy  groups  and 
l600  mesh  points,  but  these  limits  may  be  changed  very 
simply.  Since  there  are  many  operations  which  are  per¬ 
formed  repeatedly  in  the  program,  maximum  use  has  been 
made  of  subroutines. 

All  required  input  data  is  read  in  at  one  time;  both 
comment  cards  and  the  sample  problem  indicate  the  data 
needed  and  the  required  format.  Figure  l4  gives  a  flow 
chart  of  the  main  program.  A  sample  input  data,  a  glos¬ 
sary  of  variables,  and  a  sample  output  are  included  at 
the  end  of  this  appendix. 

Subroutines 

The  following  subroutines  are  called  for  in  the  pro¬ 
gram.  Other  subroutines,  such  as  computation  of  parameters 
allowing  for  anisotropic  scatter,  may  be  added  as  desired. 

Subroutine  C0N2 .  This  subroutine  computes  the  con¬ 
stants  used  in  Eqs  (16)  and  (18),  If  one  defines  the 
following  variables 
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ncos 

TSIN  %/Aj=A2 
iLr-DB^l/z  HA3 

then  Bqs  (16)  and  (18)  become 


(A ^ +A2+A^ ) PBY (1+ 1 , J )  +  fA^ +A2+A^ ) ?BY (I , J ) + (-2A2 ) PBX (I , J ) 


=r  SX(I,J)  (33) 


and 


(-A^  +A2+A^  )PBY  (I ,  J  )  -4^  )  PBY(I+ 1 ,  J  )  +  (2A2  )  PBX  (I ,  J+ 1  ) 


SX(I,J)  (34) 


where  the  variables  are  shown  in  Figure  15 • 


the  angular  position  as  shown  in  Pigura  16. 
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The  order  in  which  the  fluxes  are  calculated  is: 

L  =  9,  10,  1!,  12,  3,  4,  8,  7,  6,  5,  2,  t. 

Subroutine  FIS .  This  subroutine  determines  whether 
to  use  the  fission  cross-sections  for  fixed  fuel  or  rod 
fuel  in  computing  the  fission  source  term  SUM  for  each 
cell  in  the  core. 

Subroutine  SCAT.  To  avoid  the  necessity  of  defining 
the  scattering  cross-section  at  each  point  in  the  grid 
and  of  storing  all  this  data  in  memory,  this  subroutine 
uses  the  indices  I  and  J  of  the  point  to  locate  the  type 
of  material  associated  with  the  point.  The  appropriate 
cross-sections  are  then  stored  in  the  dummy  matrix  SCO 
to  be  used  in  the  main  program. 

Subroutine  TOTAL .  This  subroutine  duplicates  the 
function  of  subroutine  SCAT,  except  it  stores  the  total 
cross-section  for  the  point,  for  the  particular  energy 
group  the  main  program  is  using,  in  the  location  TOT. 
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Subroutine  CKZ.  Xn  ti^la  subroutine  the  an^lar 
riuxes  are  set  tc  zero  if  negative  and  are  accuaulated 
(integrated)  in  the  saatrix  p4. 

Subroutine  ROD.  In  order  to  drive  the  reactor 
critical,  soae  paraaeter  nust  be  variable  in  the  problea. 
This  V8S  accomplished  by  changing  the  volume  fractions 
of  fuel  and  cadmium  in  the  rod  fuel  regions.  This  sub~ 
routine  increases  the  amount  of  poison  if  the  reactor  is 
supercritical  vice  versa.  The  efficacy 

of  this  technique  is  doubtful,  since  the  thermal  flux 
was  greatly  depressed  in  these  regions  but  the  reactor 
remained  supercritical . 

Subroutines  INPUT?  and  INPUTS .  Those  two  subrou¬ 
tines  allow  one  to  read  and  write  binary  information  on 
tapes  7  and  8  respectively.  This  is  more  efficient 
than  using  decimal  numbers,  thus  decreasing  computational 


time . 
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THiS  PROGRAM  REQUIRES  SUBROUTINES  T0TAL»CKZ*FIS»SCAT*C0N2t IN0UT7* 
ROD*  AND  INCUTS 

2-D*  S4  CALCULATION  USING  CARLSON *S  MU{M>-BARS*  FOR  40  X  40  MESH 
DIMENSION  FS(2*3) *FCS{3) .Stl3*6) .ST(13*3J  *Pt40»40*3I » 

1SUM(8*16) *V( 3) »SX{40*40.3I *SCC(5) .PBX (40 *41 ) *PBY( 41*40  5  fKIB) *P4(40 
2*40) >P5(40*405  *DX{40) ♦DY{40)iFl (40) »F2 { 40 ) *F3 ( 40 ) ♦F4 ( 40 5  tFS ( 40 ) *F6 
3(40) ♦BXCXI 12} 

DIMENSION  QC(6) tOFtS) ,QF1 ( 3 ) »QF2 ( 3 5 ♦QC2{ 3 ) 

READ  CROSS-SECTIONS*  INCREMENTS*  VALUES  OF  NU 
INCREMENTS  IN  THE  X-DIRECTICN 
READ(5*105} (DX{ I ) *1=1.40) 

INCREMENTS  IN  THE  Y-DIRECTION 

EVERY  OTHER  INCREMENT  IS  READ,  EQUAL  TO  THE 

QUANTITY  {DY{ I )+DY( I  +  l  'D  .these  ARE  THEN  DIVIDED  BY 

TWO  AND  STORED  APPROPRIATELY.  IF  UNEQUAL  SPACING 

IS  DESIRED*  THIS  MAY  BE  CHANGED  TO  READ  EACH 

OF  THE  40  INCREMENTS. 

READ{5»105)(0Y(n  *1=1,39,2) 

DO  6000  1=2*40*2 
DY( I-1)=0Y{ I-l )/2. 

6000  DY{ I )=0Y{ I-l) 

FISSION  CROSS-SECTIONS 

FS(1*K)  CONTAINS  ROD  FUEL  FISSION  CROSS-SECTIONS. 

FS{2*K)  CONTAINS  FIXED  FUEL  CROSS-SECTIONS. 

READ(  5,106)  (  (FS(M,1C)  ,X=1*3)  *M=1*2) 

VALUES  OF  NU 
READ(5*106)(V(I),I=1,3) 

VALUES  OF  CHI 
READ(5,106) (X{ 1)91=1,3) 

CONVERGENCE  CRITERION  (EPSILON). 

READ! 5 *106) AAA 

GROUP  TO  GROUP  SCATTERING  CROSS-SECTIONS 
READ( 5*105) ( (S(M,N) ,N=1*6) *M=1,13) 

TOTAL  CROSS-SECTIONS 

READ( 5*106 ) ( (ST(M,K) ,K=1,3) ,M=1*13) 

INITIALIZE  THE  FLUX 

L  IS  THE  NUMBER  OF  THE  UNIT  ON  -WHICH  THE 
INITIAL  FLUX  IS  WRITTEN. 

READ( 5*99)L 
NUM=0 

READ{5,120) (BXCXI I )♦ 1=1,12) 

WRITE!  6, 119)  (BXCX(  I  )  ,1  =  1,12.5 

READ  {L)II*IJ,IK,({(P{I,J,K),I=1,II),J=1,IJ),<=1,IK) 

REWIND  L 

DATA  NSC*VU,VC/0, 1.  ,0./ 

D82  IS  THE  PSEUDO-ABSORPTION  CROSS-SECTION  TO 
ALLOW  FOR  VERTICAL  LEAKAGE. 

DATA  NTP,NN,NT*KTR,PI ,DB2/4#0,. 0174 5329,. 0047752/ 

DO  5005  N=l,6 
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0C{N}=S(5*N) 

5005  OF(N}=S(9»N) 

DO  5006  K=1.3 
0F1(K)=FS(1»K> 

QF2{lC)=sST{9,K) 

5006  OC2(K}-ST(5fK) 

SET  THE  THERMAL  FLUX  IN  THE  CENTER  OF  THE  CADMIUM 
CURTAIN  TO  ZERO* 

DO  1234  J=l»19 
P(12*J»3)=*0 
1234  P(13*J*3)“.0 
C  CALCULATE  SOURCE  TERMS 

5  SC£=0, 

NTP=NTP+1 

NN«NN+1 

C  CALCULATE  THE  TOTAL  FISSI'  *  SOURCE 

DO  2  Js=l*16 
DO  2  I=l»8 

CALL  FIS( I*J*FS*FCS) 

SUM{ I»J)=0. 

DO  3  K«=1.3 

3  SUM(  I»J)=SUM(  I*J)4-V{K)*FCSIIC)*P{  I  »J»K) 

2  SCE=SCE+SUM{ I*J)»DX(I)»DY( J) 

IF(KTR)35»35*6 

C  -NORMALIZE  THE  FISSION  SOURCE 

35  DO  4  I=l»8 

DO  4  J=l»16 

4  SUM(IfJ)»SUM(I.J)/SCE 

C  NORMALIZE  THE  FLUX* 

671  DO  950  K=l*3 

DO  950  1*1*40 

DO  950  J=l»40 

950  P(I*J»K)*P(I*J,K)/SCE 

C  STORE  THE  EIGENVALUE  LAMBDA 

S1*SCE 

C  CALCULATE  THE  SCATTERING-  SOURCE 

DO  40  J«l*40 
DO  40  1=1*40 
CALL  SCAT(S*SCC*I*J) 

DO  40  K*l»3 
SXU  *J*K)=0* 

DO  40  KU^lfK 
N*KU+K~( 1/KU) 

40  SX(I»J*K}=P(I*J*KU)<^SCC(N)+SX(I»J*K) 

C  CALCULATE  THE  TOTAL  SOURCE  IN  THE  CORE  FOR  EACH  CELL. 
DO  41  K=l*3 
DO  41  J*1.16 
DO  41  1=1*8 

41  SX(I»J»K)=SX(I*J*K)+X(K)*SUM( I,J) 
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C  BEGIN  INNER  ITERATIONS 

DO  32  K=l*3 

C  SET  BOUNDARY  CONDITIONS  AT  ZERO 

DO  7  1=1, AO 
P8X( I ,41 )=0. 

7  PBY{41,I}=0, 

C  L  IS  THE  INDEX  FOR  THE  AZIMUTHAL  ANGLE  PHI 

L=9 

DO  45  JD=1,40 
J=41-JD 
DO  45iD=l,40 
I=41-ID 

C  THE  TOTAL  CROSS-SECTIONS  FOR  THE  POINT! I ♦J)  ARE  DENOTED  BY  TOTIK) 
CALL  TOTAL! I ♦J,K,ST,TOT) 

CALL  CON2!A,B,C,L,DX! I ) ,OY ! J ) ,DB2 , TOT ) 

PBY! I »J}=SX! I ,J»K5-A»PBY! I+l, J) 

PBY! I,J)=!PBY! I , J ) -C*P8X ! I , J+1 ) ) /B 
P4! I ,J)=!PBY! I »J)+PBY! I+i,J) )/2. 

45  PBX! I,J)=2,*P4!I,J)-PBX! I,J+1) 

C  SET  ANY  NEGATIVE  FLUX  TO  ZERO. 

DO  1000  I=;l,40 
DO  1000  J=l,40 
IF!P4! I , J) ) 1001,1000,1000 
1001  P4!I»J)=0. 

1000  CON'I  INUE 

THE  MATRICES  F  ARE  USED  TO  STORE  THOSE  ANGULAR 
FLUXES  reflected  ACROSS  A  LINE  OF  SYMMETRY  INTO 
ANOTHER  ANGLE. 

DO  49  1=1,40 
F3!I)=P8Y!1,I) 

49  F2! I )=PBX! I  ,1) 

PBX!I,1)  FOR  L=9  IS  P8X!I,1)  FOR  L=8 
P8Y!1,J)  FOR  L=9  IS  PBY!l,J)  FOR  L=i2 
L  =  10 

DO  51  JD=1,40 
J=41-JD 
DO  51  10=1,40 
I=41-ID 

CALL  TOTAL! I ,J,K,ST,TOT) 

CALL  GON2(A,B,C»L,DX!n ,DY ! J )  ,DB2 ,TOT ) 

PBY! I ,J)=SX! I ♦J,K)-A*PBY! I+1,J) 

PBY! I,J)=!PBY! I , J ) ~C*P8X ! I , J+1 ) ) /B 
P5! I ,J)=!PBY! I ,J)+PBY! I+1,J) )/2. 

51  PBX!  I  ,J)=2.^^P5!  I  ,J)-PBX!  I,J+1) 

CALL  CKZ!P4,P5) 

C  PBX!I,1)  FOR  L=10  IS  PBX!I,1)  FOR  L=7 
DO  55  1=1,40 
55  FI! I )=PBX! 1 ,1) 
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00  219  J0=lt40 

J=41-JD 

DO  219  1=1*40 

CALL  TOTAL! I »J»K»ST,TOT) 

CALL  C0N2(A*B.C*L»DX( I ) tOY { J 1 *082 ♦ TOT > 

PBY(  I+l  •J)  =  (SX{  I  »J*K)-C»PBX{  I  tJ+D-B^POY!  I  *J)  )/A 
P5{ I tJ)=(PBY( I+l .J)+PBY( I tJ) )/2* 

219  PBX{  I»J)=2.»P5{I*J3-PBX(  r^,J+l) 

CALL  CKZ!P4,P5) 

C  PBXCUl)  FOR  L=ll  IS  PBXd,!)  FOR  L  =  6 

DO  122  I=l»40 
122  F4tI)=PBX(I*l) 

L=12 

DO  57  J=l*40 
57  PBY! l.J)=F3{ J) 

DO  59  JD=1.40 

J=41-JD 

DO  59  1=1*40 

CALL  TOTAL! I ♦J»K.ST»TOT) 

CALL  CON2(A.8*C.L*DX! I ) ,DY1J) *082* TOT) 

P8Y! I+1*J)=!SX! I *JtK)~C*P8X! I » J+1 ) -B+P8Y ! I .J) )/A 
P5! I ♦J)=!P8Y! I+1*J)+PBY! I *J) )/2. 

59  PBX! I ♦J)=2.*P5! I fJl-PBX! I ,J+1) 

CALL  CKZ!P4*P5) 

C  PBX! 1,1)  FOR  L=12  IS  PBX!I,1)  FOR  L=5 
DO  64  1=1,40 
64  F3!  n=P8X!  I  ,1 ) 

L=3 

DO  73  1=1,40 
PBX! I,41)=0, 

73  PBY!41,I)=0, 

DO  10  JD=1,40 
J=41-JD 
DO  10  10=1,40 
1=41-10 

CALL  TOTAL! I ,J,K,ST,TOT) 

CALL  C0N2(A,B,C,L,DX( I ) ,DY ! J ) ,032 ,TOT ) 

PBY! I ,J)=!SX! I ,J,K>-A*PBY! I+1,J)-C*PBX! I ,J+1) )/B 
P5! I ,J)=! PBY! I ,J)+PBY! I+1,J) )/2. 

10  PBX!I,J)=2*»P5!I,J)-PBX!I,J+1) 

CALL  CKZ!P4,P5) 

DO  13  1=1,40 
13  F5!I)=PBX!I,1) 

L=4 

DO  19  JD=1,40 

J=41-JD 

DO  19  1=1,40 

CALL  TOTAL! I ,J,K,ST,TOT) 

CALL  C0N2!A,B,C>L,DX! I ) ,DY ! J ) ,DB2 ,T0T ) 
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PBY( I+1*J)={SX( I .J,K)-C*PBX( I .J+1)~B*PBY{ I tJ) )/A 
P5{ I *J)={PBY{ I+l, J)+PBY{ I ♦J) )/2# 

19  PBX(  I*J)=2**P5{I.J)-'PBX{  I»J+1) 

CALL  CKZ(PA»P5) 

DO  22  1=1.40 
22  F6( I )=PBX( I.l) 

L=8 

DO  67  1=1.40 
PBY{41. I )=0. 

67  PBX( I .1)=F2( I ) 

DO  126  J=1.40 
DO  126  ID=1»40 
I=41-ID 

CALL  TOTAL( I .J.K.ST.TOT) 

CALL  C0N2(A.B.C.L.DX{ I ) .DY ( J ) .082 .TOT ) 

PBY{ I .J)=(SX( 1 .J.K)-C*PBX{ I .J)-A»PBY( I+l.J) )/B 
P5(I»J)=(PBY( I+1.J)+PBY( I.J) >/2. 

126  PBXt  I  .J+1)==2.*P5(  I .  J)-PBX(  I  .J) 

CALL  CKZ(P4.P5) 

C  PSY(l.I)  FOR  L=8  IS  PBYtl.I)  FOR  L=5 

DO  69  1=1.40 

69  F2( I )=PBY( 1. I ) 

L  =  7 

DO  70  1=1.40 

70  PBX{I.1)=F1(I) 

DO  225  J=l»40 
DO  225  10=1.40 
I=41-ID 

CALL  TOTAL! I .J.K.ST.TOT) 

CALL  C0N2{A.B.C.L.DX( I ) .DY ( J ) ,DB2 . TOT ) 

PBY( I .J)  =  (SX( I .J.K)-C-»PBX{ I  .J)~A»PBY( I  +  l.J) )/B 
P5( I .J)=(PBY( I+l.J )+PBY( I.J) )/2. 

225  P8X{  I  .J+1)=2.^^P5(  I  .J)-PBX(  I  .J) 

CALL  CKZ(P4.P5) 

L=6 

C  PBYd.J)  FOR  L=7  IS  PBYd.J)  FOR  L=6 

DO  71  J=1.40 

71  PBX( J.l )=F4( J) 

DO  127  J=1.40 
DO  127  1=1.40 

CALL  TOTAL! I .J.K.ST.TOT) 

CALL  CON2!A.B.C»L.DXd ) .DY! J) .DB2.TOT) 

PBY! I+l ♦J)=SX! I .J.K)-C*PBX! I ♦ J) 

PBY!  I  +  l.J)  =  !PBY!  I+l.J) -B^^PBY!  I  .J)  )/A 
P5! I .J)=!PBY! I ♦J)+PBY! I+l.J) )/2. 

127  PBX! I ♦J+1)=2.*P5! I .J)-P3X! I.J) 

CALL  CKZ!P4.P5) 

L  =  5 

DO  72  1=1.40 
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PBY( !♦ I )=F2{ I ) 

72  PBX{ 1 ,1 )  =  F3(  I  ) 

DO  227  J=i»40 
DO  227  I=l»40 
CALL  TOTAL{ I ,J»K,ST»TOT) 

CALL  CON2(A*8*C*L,DX( I ) *DY ( J ) *DB2 *TOT ) 

PBY( I+1»J)=SX( I»JfK)~C*PBX{I«J) 

PBY{ I+1»J)-(PBY( I+1,J)-B*PBY{ I*J) )/A 
P5(I »J)={PBY{ I *J)+PBY( I+1*J) )/2* 

227  PBX( I ♦J+1)=2.*P5( I*J)-PBX( I .J) 

CALL  CKZ{P4,P5) 

L=2 

DO  23  1=1*40 
P8Y(41,n=0. 

23  PBX( !»1)=F5(  I  ) 

DO  25  J=1.40 
DO  25  10=1*40 
I=41-ID 

CALL  TOTAL( I »J.K»ST.TOT) 

CALL  CON2(A,B*C*L*DX{ 1 ) tDY ( J ) *DB2 *101 ) 

PBY( I .J)=(SX( I ♦J*K)-C*PBX( I fJI-A^PBYC I+1*J) ) /B 
P5( i »J)=(PBY( I+1»J)+PBY( I tJ) )/2. 

25  PBX( I*J+1)=2.*P5{ I*J)-PBX{ I*J) 

CALL  CKZ(P4,P5) 

L=1 

DO  66  1=1*40 
66  PBX(  I*l)=F6n  ) 

DO  27  J=l*40 
DO  27  1=1*40 

CALL  TOTAL( I *J*K*ST*TOT) 

CALL  C0N2(A*B»C*L*DX( I ) »DY ( J ) *DB2 *TOT ) 

PBY( I+l ,J)=(SX( I *J*K)-C*PBX{ I ♦J)-B*PBY( I *J) ) /A 
P5( I *J)={PBY( I *J)+PBY( I+1,J) )/2. 

27  PBX( I*J+1)=2**P5( I *J)-P8X( I *J) 

C  CALCULATE  TOTAL  FLUX  AT  EVERY  POINT 
DO  1006  J=l*40 
DO  1006  1=1*40 
IF(P5{ I ,J) ) 1007,1006,1006 
1007  P5( 1 *J) =0* 

1006  P4( I ,J) ={P4( I ,J)+P5( I*J) )/12. 

IF(K-3> 1235*1236,1236 

1236  DO  1237  J=l*19. 

DO  1237  1=12,13 

1237  P4( I  ,J)  =  .0 

C  IMPROVE  SCATTERING  SOURCE  FOR  NEXT  GROUP 

1235  DO  32  J=l*40 
DO  32  1=1*40 

C  STORE  THE  APPROPRIATE  CROSS-SECTIONS  IN  SCC. 
CALL  SCAT(S,SCC*I * J) 
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DO  33  L=K»3 
N=K+L-( 1/K) 

33  SX( I ♦J»L)=SX( I #J*L)+SCC(N)*(P4{ I , J ) -P ( I , J ,K ) ) 

32  P{I»J»K)=P4( I»J) 

C  END  OF  LOOP 

C  WRITE  THE  TOTAL  FLUX  ON  TAPES  7  AND  8  ALTERNATELY. 

C  THIS  ALLOWS  RESTARTING  THE  PROGRAM  WITH  THE  BEST 

C  AVAILABLE  FLUX  VALUES.  TWO  TAPES  ARE  USED  IN  CASE 

C  THE  PROBLEM  TERMINATES  BEFORE  THE  LAST  WRITE 

C  COMMAND  CAN  BE  EXECUTED. 

IF(NUM) 1162*1162.1163 

1162  NQQ=8 
NUM=1 

GO  TO  1164 

1163  NQQ=7 
NUM=0 

1164  11=40 
IJ=40 
IK=3 

WRITE(NOQ) II»IJ»IK»(((P(I.J.K).I=1.II)»J=1»IJ)*K=1,IK) 

END  FILE  NQQ 

REWIND  NQQ 

WRITE! 6 *109) NQQ 

KTR=1 

C  CALCULATE  IMPROVED  FISSION  SOURCE 

GO  TO  5 

C  HAVE  INNER  ITERATIONS  CONVERGED 
6  R«(SCE/S1)-1. 

I F{ABS(R) -AAA >34.34.357 

C  WRITE  THE  FLUX  EVERY  25  ITERATIONS  WHETHER  OR  NOT 
C  CONVERGENCE  HAS  BEEN  ACHIEVED. 

357  IF{NTP-25)35.358.358 

358  WRITE(6,121)R 
WRITE(6.100)SCE 
WRITE(6.107)NN 
WRITE{6,98) 

WRITE(6.102) 

WRITE(6.103) ( { I . J. (P( I .J.K) .K=1.3) . 1=1.40) .J=1.40) 

NTP=0 

GO  TO  35 

34  R=SCE-1. 

C  HAVE  OUTER  ITERATIONS  CONVERGED 
IFIABSIR)-. 01)36. 36.37 
36  WRITE(6.104) 

NT=0 

A=P(1.1.3) 

DO  970  K=1.3 
DO  970  J=l*40 
DO  970  1=1.40 


-m 


on  on 


970 

P!  I 

♦  J 

.K) 

=  P!  I 

.J,.<)  /A 

KTR 

=“ 

1 

37 

IF! 

NT 

)60 

,377 

.378 

378 

IF! 

NT 

-50 

)61 . 

377,377 

377 

WRI 

TE 

!  6 . 

101) 

DB2 

WRI 

TE 

!  6 . 

100) 

SCE 

WRI 

!E 

!6. 

107) 

NN 

WRI 

TE 

!6. 

98) 

WRI 

TE 

!6« 

102) 

WRITE(6*103) { ( I » J» (P( I tJ»K) .K=l»3) ♦ I=l»40) ♦J=l*40) 
NT=0 

IF(DB2)81»65»81 
81  IF{KTR)80*65»61 

61  IF{NSC)5000. 5000. 5001 

5000  IF(R)5.377,5001 

ADJUST  THE  CONTROL  ROD  BY  VARYING  THE  VOLUME 
FRACTIONS  OF  FUEL  AND  CADMIUM. 

5001  CALL  R0D(FS.S.ST.VU9,VC.NSC.R.QC.QF,QF1,QF2»0C2) 
WRITE(6,108) VU 

IF(VU-. 05)60.5011.5011 
5011  KTR=0 

NT=NT+1 
GO  TO  5 

AFTER  CONVERGENCE.  CALCULATE  THE  FLUX  WITH  AN 
INFINITE  HEIGHT  CORE. 


80 

65 

98 

99 

100 
101 
102 

103 

104 

105 

106 

107 

108 
109 
113 

119 

120 
121 
60 


C 


DB2=0. 

GO  TO  35 
WRITE(6.113) 

FORMAT! IH  ) 

FORMAT! 18) 

FORMAT! IH  ♦21HTOTAL  FISSION  SOURCE=  .£9.3 ) 

FORMAT! 1H1.14HBUCKLING  TERM=.E10.3) 

FORMAT !39H  X  Y  GROUP  1  GROUP  2  GROUP  3) 
F0RMAT!1H  .2  I  3 . IX . 3E1 1 . 3 ) 

FORMAT! 1H1.28HCONVERGED  FLUX  VALUES  FOLLOW) 
FORMAT!6E12.6) 

FORMAT!3E12.6) 

FORMAT! IH  .21HNUMBER  OF  I TERAT IONS=  .  1 8  ) 

FORMAT! IH  .24HVOLUME  FRACTION  OF  FUEL=.E10.3) 
FORMAT! 17H  FLUX  IS  ON  TAPE  .13) 

FORMAT! IH  .ISHNORMAL  TERMINATION) 

FORMAT! 22A6) 

FORMAT! 12A6) 

FORMAT! IH  .12H!SCE/S1 )-!.=. F15. 5) 

STOP 

END 

SUBROUTINE  CON2! A.&.C.L.DX.DY.DB2.TOT) 

SEPARATE  LOOP  TO  CALCULATE  CONSTANT  COEFFICIENTS 
XL=L 
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IF(L-4) 14*14tl414  * 

14  A=2*»XL-1. 

A=A*3. 14159/4, 

C  USE  THIS  VALUE  FOR  L=l»4 

ETA=. 8703883**2 
ETA=1,-ETA 
ETA=SORT(ETA) 

GO  TO  1415 
■1414  A=3. 14159/4, 

A=A*XL 

A=A-9. *3. 14159/8, 

C  USE  THIS  VALUE  FOR  L=5»12 

ETA=,3481553**2 
ETA=1,-ETA 
ETA=SORT(ETA) 

1415  AC=COS(A) 

AS=SIN{A) 

A1=FTA*AC/DX 

A2=bTA*AS/DY 

A3=(DB2+TOT)/2, 

IF(L-4) 1416,1416,1417 

1416  IF(L--3)30*29,29 

1417  IF(L-9)30»29»29 

29  A=A1-A2+A3 
B=-A1-A2+A3 
C=2,*A2 
RETURN 

30  A=A1+A2+A3 
B=-A1+A2+A3 
C=-2,*A2 
RETURN 

END 

SUBROUTINE  F I S ( I ♦ J, FS, FCS ) 
DIMENSION  FS{2,3) fFCSI 3) 

IF( I ,6E,3)GO  TO  3000 

IF( J,GE,13)GO  TO  3001 

IF{ J,GE,5,AND,J,LT,9)GO  TO  3001 

3000  DO  3002  KK=1,3 

3002  FCS(KK)=FS(2»KK) 

RETURN 

3001  DO  3003  KK=1»3 

3003  FCS{'<K)=FS(  IfKK) 

RETURN 

END 

SUBROUTINE  CKZ(P4,P5) 

DIMENSION  P4{40,40) »P5(40»40) 

DO  1002  I=l»40 
DO  1002  J=l,40 
IF(P5( I ,J) ) 1003»1002,1002 
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1003  P5{I»J)=0. 

1002  P4( I ,J)=P4( I ,J)+P5{ I,J) 

RETURN 
END 

SUBROUTINE  SCAT ( S*SCC» I  * J ) 

DIMENSION  S( 13*6) .SCC(6) 

THIS  SUBROUTINE  SEARCHES  THROUGH  LOGIC  STATEMENTS  TO  LOCATE  THE 
POINT! I tJ)  IN  THE  MODEL*  AND  STORES  THE  APPROPRIATE  SCATTERING 
CROSS-SECTIONS  IN  THE  MATRIX  SCC. 

IF( I.GE.37>60  TO  3015 
IF( I.GE.17.AND.J.GE*33)G0  TO  3015 
IF( I.GE.17,AN0.I.LE.18.AND.J*GE.29)G0  TO  3015 
IF( I.LE.16.AN0.J.GE*29)G0  TO  3017 
IF( I.LE*18.AND.J.GE.25*AND.J.LE.28)GO  TO  3019 
IF( I.GE.17*AND.I*LE.18*AND,J*LE.28)G0  TO  3019 
IF! I.GE.21.AND.I.LE.35.AND.J.LE*20)GO  TO  3019 
IF! I.GE.23.AND.I.LE.35.AND.J.LE.28)G0  TO  3019 
IF! I.EQ.36.AND.J*LE.32)G0  TO  3021 
IF! I.GE.19.AND.I*LE.36.AND.J.GE.29)GO  TO  3021 
IF! I*GE,19«AND.I.LE.22.AND.J«GE.21)G0  TO  3021 
IF! I.GE.19.AND.I.LE*20.AND*J.LE.20)GO  TO  3023 
IF!  I*GE.i5.AND*ULE*.16.AND.J.LE.24)G0  TO  3025 
IF! I.GE*11,AND.I.LE.14*AND*J.LE.20)GO  TO  3027 
IF! I*GE,ll*AND.I.LE.14.AND.J.LE.24)GO  TO  3029 
IF!  I.GE.,9.AND*I.LE*10.AND*J*LE.24.)GO  TO  3031 
IF! J*GE*21«AND#J.LE.24.AND.I.LE.8)GO  TO  3033 
IF! J.GE.17.AND.J*LE.20*AND.I*LE.8)GO  TO  3035 
IF!I.GE,3.AN0,ULE.8,AND,J*LE.16)G0  TO  3037 
IF!J.LE.4)G0  TO  3037 
IF! J.LE*12oAND.J*GE.9)GO  TO  3037 
GO  TO  3039 
C  CONCRETE 

3015  DO  3016  N=l*6 

3016  SCC!N)=S! 13*N) 

RETURN 

C  WATER 

3017  DO  3018  N=l*6 

3018  SCC!N)=S! 12*N) 

RETURN 

C  AIR 

3019  DO  3020  N=l»6 

3020  SCC!N)=S! 11*N) 

RETURN 

C  TEST  CELL  WALL 

3021  DO  3022  N=l*6 

3022  SCC!N)=S! 10*N) 

.  RETURN 

C  AL-H20  4 

3023  DO  3024  N=l*6 


51 


r>  r>.  o 


GNE/PHYS  65-10 


302+  SCC(N)=S(4,N) 

RETURN 

C  AL-H20  3 

3025  DO  3026  N=l*6 

3026  SCC{N)=S(3*N) 

RETURN 

C  CADMIUM 

3027  DO  3028  N=l*6 

3028  SCC(N)=S(5.N) 

RETURN 

C  ALUMINUM 

3029  DO  3030  N=1.6 

3030  SCC(N)=S(6*N) 

RETURN 

C  AL-H20  1 

3031  DO  3032  N=l»6 

3032  SCC(N)=S( 1*N) 

RETURN 

C  AL-H20  2 

3033  DO  3034  N=lf6 

3034  SCC{N)=S(2»N) 

RETURN 

C  8E-H20  REFLECTOR 

3035  DO  3036  N=l*6 

3036  SCC{N)=S(7fN) 

RETURN 

C  CORE 

3037  DO  3038  N=i*6 

3038  SCCIN)=S(8.N) 

RETURN 

C  ROD  FUEL 

3039  DO  3040  N=l»6 

3040  5CC(N)=S(9»N) 

RETURN 

END 

SUBROUTINE  TOT AL ( I  * J*K ♦ ST ♦ TOT ) 

DIMENSION  ST(13»3) 

THIS  SUBROUTINE  ALSO  SEARCHES  THROUGH  LOGIC  STATEMENTS  TO  LOCATE 
THE  POINT! I *0)*  AND  STORES  THE  APPROPRIATE  TOTAL  CROSS-SECTION 


FOR 

ENERGY 

GROUP 

!K) 

IN  THE 

LOCATI 

ON 

TOT 

IF(  I 

•  GE  • 

37 

)GO 

TO 

401 

5 

IF(  I 

•  GE  • 

17 

.AND 

•  J  i 

.GE. 

33)GO 

TO 

4015 

IF(  i 

•  GE  • 

17 

•  AND 

•  u 

,LE. 

18. AND 

•  J 

.GE.29 

)GO 

TO 

4015 

IF{  I 

•  LE  • 

16 

.AND 

•  J  i 

.GE. 

29)  GO 

TO 

4017 

IF{  I 

.LE. 

18 

.AND 

•  J  i 

.GE. 

25. AND 

•  J 

.LE.28 

)GO 

TO 

4019 

IF(  I 

•  G£  • 

17 

•  AND 

.LE. 

18. AND 

m  J 

.LE.28 

>GO 

TO 

4019 

IF(  I 

•  GE  • 

21 

•  AND 

.1. 

.  le  . 

35. AND 

.  Lw  .  20 

)GO 

TO 

4019 

IF!  I 

.GE. 

23 

.AND 

•  I . 

.LE. 

3 5. AND 

•  J 

.LE.28 

)GO 

TO 

4019 

IF!  I 

.EO. 

36 

.AND 

•  J  1 

.LE. 

32)GO 

TO 

4021 
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IF( I.GE.19.AND.I.LE.36*AND*J*GE*29)G0  TO  4021 
IF(I.GE.19.AND.r.LE.22.AND*J.GE.21)GO  TO  4021 
rF( I,GE,19.ANO.T*LE.20*AND*J*LE.20)GO  to  4023 
IF(  I.GE,r5.AND.I.LE«16.AND.J.LE.24.)60  TO  4025 
IF(I.6E.11.ANO.I.LE.14.ANO*J.LE.20)GO  TO  4027 
IF{I.GE*11.ANO.I.LE.14.AND.J.LE.24)GO  TO  4029 
IF(I*GE,9.AND.I.LE.10.AND*J.LE*24)GO  TO  4031 
IF( J.GE*21.AND.J*LE*24.AND.I.LE.8)G0  TO  4033 
IF( J.GE,l7,AND*J.LE*20.ANO.i.LE.8)GO  TO  4035 
IF{ I.GE*3.AND.I.LE.8.AND.J.LE.16)GO  TO  4037 
IF{J,LE.4)G0  TO  4037 
IF( J.LE.12.AND.J.GE*9)GO  TO  4037 
4039  TOT=ST(9.K) 

RETURN 

4015  TOT=ST(13*K) 

RETURN 

4017  TOT=ST(12*K) 

RETURN 

4019  TOT=ST(ll*K) 

RETURN 

4021  T0T=ST(10*K) 

RETURN 

4023  T0t=ST{4,K) 

RETURN 

4025  TOT=ST(.3,*K) 

RETURN 

4027  TOTsSTJ5>K) 

RETURN 

4029  TOT«ST(6*K) 

RETURN 

4031  TOT=ST(l»K) 

RETURN 

4033  TOT=ST(2*K) 

RETURN 

4035  TOT=ST(7fK) 

RETURN 

4037  fOT»ST{8*K) 

RETURN 

END 

SUBROUTINE  ROD( FS*S *STt VUt VC*NSC»R»QC»QF *QF1 »QF2 *QC2 ) 
DIMENSIONFS(2t3) *S( 13f6) *ST( 13*3 ) *QC (6) *QF( 6 ) »QF1 ( 3 ) ♦QF2 ( 3 ) ♦ 
1QC2(3) 

NSC=NSC+1 

IF(R)1000.1001»1001 

1000  0=1. 

GO  TO  1002 

1001  Q=-l. 

1002  VU*{  .S^^NSO^Q+VU 

vc*i.-vu‘ 
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DO  1003  •N=l*6 

0F(N)  CONTAINS  THE  ORIGINAL  SCATTERING  CROSS-SECTIONS 
FOR  THE  ROD  FUEL.  AND  QC(N)  THOSE  FOR  CADMIUM. 

1003  S{9»N)=QF(N)*VU+QCrN)*VC 
DO  1004  K=1.3 

FSd.K)  IS  ROD  FUEL.  FS(2.K)  IS  FIXED  FUEL. 

THE  INITIAL  FISSION  CROSS-SECTIONS  ARE  STORED  IN  QFl. 
FS( l.K)=VU*QFl (K) 

THE  INITIAL  TOTAL  CROSS-SECTIONS  FOR  ROD  FUEL  AND 
CADMIUM  ARE  STORED  IN  OF2  AND  0C2  RESPECTIVELY. 

1004  ST(9»K)=OF2(K)*VU+QC2{K)*VC 
RETURN 
END 

SIBMAP  INOUT7  4 

ENTRY  .UN07. 

.UN07.  PZE  UNIT07 

UNIT07  FILE  .A( 1 ) .BIN. INOUT.BLK=256.NOL 1ST. READY 
END 

SIBMAP  INOUT8  4 

ENTRY  .UN08. 

.UN08.  PZE  UNIT08 

UNIT08  FILE  .C(  1 )  .BIN.  INOUT.BLIC==256.NOLIST. READY 
END 
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C  SAMPLE 

INPUT  DATA 

FOR  S4C40 

2.425 

2.425 

2.57 

2.57 

2.57 

2.57 

2.57 

2.57 

.2285 

.2285 

.051 

.051 

.051 

.051 

.96 

.96 

6.985 

6.985 

.6905 

.6985 

5.65 

5.65 

22.5 

22.5 

22.5 

22.5 

22.5 

22.5 

22.5 

22.5 

22.5 

22^5 

22.5 

22.5 

22.5 

12.7 

22.5 

22.5 

22.5 

22.5 

.23127E+01 

.23127E+01 

.38545E+01 

i38545E+0l 

.38545E+01 

.38545E4-01 

.38545E+01 

.38545E+01 

.38545E+01 

.38545E+01 

.8255E-00 

.8255E-00 

.69855E+01 

.69855E+01 

.635E+01 

.635E+01 

.225E+02 

.225E4-02 

.225E+02 

.225L+02 

.15557E-03 

.27244E-02 

.446E-61 

.904E-04 

.127E-02 

.304E-01 

.261E+01 

.243E+01 

.246E+01 

.lE+01 

.37E-05 

.OE~00 

.001 

.222E-00 

.245E-02 

.156E-03. 

.119E-00 

.422E-02 

.1.765E-00 

.2759E-00 

.186E-01 

.135E-05 

.373E-00 

.320E-01 

•878E-00 

.249E-0b 

.108E-0i 

.775E-06 

.250E-00 

.186E-01 

.539E-00 

.227E-00 

.524E-02 

.363E-06 

.151E-00 

.372E-02 

•255E-00 

.277E-00 

.2747E-03 

O.E-00 

.226E-00 

.462E-03 

•322E-00 

.2127E“-00 

.3493E-03 

O.E-00 

.856E-01 

.5989E-03 

.84E-01 

.446E-00 

.122E-01 

.266E-06 

.720E-00 

.194E-01 

.930E-00 

.378E+00 

.484E-01 

.3571E-05 

.8439E+00 

.8331E-01 

.217E+01 

•38652E-00 

.50856E-01 

.3751E-05 

.888E-00 

.88886E-01 

.2297E+01 

.274E-00 

.294E-02 

.173E-06 

.571E-00 

.243E-02 

.667E~00 

.92E-04 

.589Er-09 

•OE-00 

.249E-03 

.0E~00 

.381E-03 

.483E-00 

.791E-01 

.585E-05 

.132E+61 

.136E-'00 

.351E+01 

•243E-00 

.197E-01 

.138E-05 

.442E-*00 

.3E-01 

.828E-00 

.22452E-00 

.12426E-00 

.1891E-00 

.29478E-00 

.40551E-00 

.892E-00 

.26036E-00 

.26945E-00 

.542E-00 

.232457E-00 

.15555E-00 

.268E-00 

.2868E-00 

.4285E-00 

.1355E+03 

.2134E~00 

.8687E-01 

.962E-01 

.46067E-00 

.73928E-00 

.931E-00 

.427E-00 

.93337E-00 

.871E-01 

.43783E-00 

.98225E-00 

.23746E+01 

.277E-00 

.57637E-00 

.787E-00 

4 

.9525E-04 

.2492E-03 

.438E-03 

.5627E-00 

.14627E+01 

.353E+01 

.2634E-00 

.4734E-00 

.848E-00  • 

8 

ALL  INPUT  DATA  HAS  BEEN  READ  CORRECTLY 
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Gloa»ary  of  Coaputor  Prograa  Sybols 


A 

• 

• 

AziHuthal  aniflo  |^| 

AAA 

: 

Convergenea  critorlon | £ | 

DB2 

X 

Vortical  laakafo  tor*{  DBg^) 

DX 

t 

Incranent  in  tha  x-dlractlon 

DY 

• 

• 

Incramant  in  tha  y-diraetlon 

FS 

• 

• 

Flsalon  croas-aactlon 

NN 

• 

• 

Numbar  of  iterations 

P 

• 

• 

Total  neutron  flux  [N{r)] 

PBX 

• 

• 

Ancrular  flux  (N^.j) 

PBY 

• 

• 

Angular  flux|,fsjj  j| 

S 

: 

Scattering  cross-section 

SCE 

• 

• 

Total  fission  source  (^| 

ST 

• 

• 

Total  cross-section 

r 

SUM 

• 

• 

Fission  source  in  one  cell 

sx 

« 

• 

Scattering  source  in  one  cel] 

L 

V 

• 

• 

Average  nusiber  of  neutrons  per 

fission 

VC 

• 

• 

VoluBie  fraction  of  cadmiim  in 

rod  fuel 

vu 

• 

• 

Volume  fraction  of  fuel  in  rod 

fuel 

X 

• 

• 

Group  source  |^| 

\v] 
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BUCKLING  TLRM=  0.478t-02 
TOTAL  FISSION  SOUR CE=0 . 1705  02 
NUMBER  OF  ITERAriONS=  5 


X 

Y 

GROUP  1 

liROUP  2 

GROUP  3 

1 

1 

.  0,688i:-02 

0.373E-03 

0.6536  00 

2 

1 

-A 

0.740£-02 

0.3986-03 

0.7206  00 

3 

1 

0.772E-02 

0.418E-03 

0.7406  00 

4 

1 

0.781E-02 

0. 4256-03 

0.7386  00 

5 

1 

0.767E-02 

0. 4166-03 

0.7376  00 

6 

1 

0.690E-02 

0.3746-03 

0.6656  00 

7 

1 

0.572C-02 

0.3116-03 

0.5466  00 

8 

1 

0.3616-02 

0.1866-03 

•  0.3996-00 

9 

1 

0.198E-02 

0.8976-04 

0. 2986-00 

10 

1 

0.1766-02 

0.857E-04 

0. 2756-00 

11 

1 

0.1646-02 

'6.8246-04 

0.3596-01 

12 

1 

0. 159E-02 

0.799E-04 

0. 

13 

1 

0.1556-^02 

0. 7766-04 

0. 

14 

1 

0.1516-02 

0.7556-04 

0. 1356-03 

15 

1 

0.1196-02 

0.6396-04 

.0.  5386-01 

16 

1 

0.7176-03 

0.4406^-04 

0. 1216-01 

17 

1 

0.5096-03 

0.323e--04 

0.4466-02 

18 

1 

0.4246-03 

0.2626-04 

0.4366=^02 

19 

1 

0.348E-03 

0.2116-04 

0. 3945-02 

20 

1 

0. 2786-03 

0. 1796-04 

0.3216-^02 

21 

1 

0.2316-03 

0.150^-04 

0. 239E-02 

22 

1 

0.200l— 03 

0. 128E-04 

0. 1926-02 

23 

1 

0.1341; -03 

0.8306-05 

0. 120E-02 

24 

1 

0.5616-04 

0.3436-05 

0. A136-03 

25 

1 

0.195,6-05 

0.3706-06 

0. 1076-04 

26 

1 

0.2376-06 

0.4076-07 

0.8226-05 

27 

1 

0.8706-05 

0.4426-06 

0.6656-04 

28 

1 

0.4606-05 

0.324E~06 

0. 1106-04 

29 

1- 

. 0.8236-06 

0.3636-07 

0. 9166-05 

30 

1 

0.3646-05 

0.2226-06 

0, 3976-04 

31 

1 

0.192E-05 

0.1486-06 

0.7866-05 

32 

1 

0.4016-06 

0.6166-07 

0.6566-05 

33 

1 

0.3876-05 

0.1956-06 

0. 2556-04 

34 

1 

0.3076-06 

0. 4276-07 

0. 7826-05 

35 

1 

0.3146-06 

0.8746-08 

0. 1596-04 

36 

1 

0.310E-05 

0.9816-07 

0. 1946-04 

37 

1 

0. 1406-06 

0.1286-07 

0. 418E-05 

38 

1 

0.7766-06 

0.5386-07 

0.9516-05 

39 

1 

0.1226-06 

0. 1106-0,7 

0. 384E-05 

40 

1 

0.2796-06 

0.2666-07 

0.  864E-05 

1 

2 

0.6546-02 

0.3546-03 

0.6226  00 

2 

2 

0.7096-02 

0.3826-03 

0.685E  00 

3 

2 

0.755U-02 

0.4096-03 

0.7226  00 

4 

2 

0.7806-02 

0.4236-03 

0.7476  00 

5 

2 

0.7666-02 

0.4166-03 

0. 7376  00 
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8  40 

0. 123L'-08 

0.416C-09 

0. 

197c -06 

9  40- 

0. 300,5-08 

0. 125E-08 

0. 

,i3ir-n4 

10  40 

0.645c5-07 

0.311L-03 

0. 

8.6  5  •: -05 

11  40 

n.l37t-06 

0.533h-0'.' 

0. 

17  7.- -05 

12  40 

0.  l/ll:-06 

0.619t-08 

0. 

i0ih*-04 

13  40 

0.194E-06 

0.696E-08 

0. 

5886-04 

14  40 

0.209E-06 

0.730E-08 

0. 

i86c-04 

15  40 

0.e96e-07 

0.368C-08 

0. 

lC5t.-()4 

16  40 

0.r39t:-07 

0.858f:-09 

0. 

1036-04 

17  40 

0.391E-08 

0.23'3E-0  ? 

0. 

1586-05 

18  40 

0.103E-07 

0.734t-09 

0. 

1446-05 

19  40 

0.105L-07 

().224C-03 

0. 

511C-05 

.20  40 

0.a89L-07 

0.157t-03 

0. 

3016-04 

21  40 

0.151E-07 

0. 129E-08 

0. 

3666-05 

22  40 

0.212C-07 

0.779C-09 

0. 

2636-05 

23  40 

0.548e-08 

0.243t-00 

0. 

8436-06 

24  40 

0.395b-08 

0.215E-09' 

0. 

326.6-06 

25  40 

0.599C-08 

0.345t-09 

0. 

8516-06 

26  40 

0.244E-08 

0.641E-09 

0. 

:i92F.-06 

27  40 

0.633C-08 

0.349L-09 

0. 

869L--06 

28  40 

0.132E-09 

0.22lfc-09 

0. 

;i826-06 

29  40 

0.101 8-07 

0.390C-09 

0. 

888L-06 

30  40 

0,213fc-09 

0.26811-09 

0. 

3776-06 

31  40 

0.164l“07 

0.504f -09 

0. 

8656-06 

32  40 

0.253E-09 

0.930t-09 

0. 

555-:-06 

33  40 

0.192L-07 

0.651C-09 

0. 

9901 -06 

34  40 

0.220L-09 

0.406C-09 

0. 

422L-06 

35  40 

0. 165L-07 

0.536U-09 

0. 

1386-05 

36  40 

0.29  IE -09 

0.2876-09 

0. 

7706-06 

37  40 

0. 115E-07 

0.5266-09 

0. 

1256-05 

38  40 

0.  142:: -09 

0.1846-09 

0, 

a53l-06 

39  40 

0.  115c-t)7 

0.662h-09 

0. 

12i'.-05 

40  40 

0. 133C-09 

0.219t-09 

0. 

4636-06 

VOLUME 

FRACTION  OF 

FiJtL=  0.500 

:  oc 

Flux  IS 

ON  TAPl  8 

FLUX  IS 

ON  TAPE  7 

FLUX  IS 

ON  FAPc  8 

. 

FLUX  IS 

ON  TAPl  Y 

FLUX  IS 

UN  TAPl  8 

VOLUME 

FKACriON  OF 

FUn-L=  0.25OE-0U 

FLUX  IS 

ON  TAPl  7 

FLUX  IS 

ON  FAPu  8 

VOLUME 

FRACTION  OF 

FUcL=  0.1256-00 
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Appendix  B 
The  GAMrl  Code 

The  GAN-I  code  (Refs  10,11)  is  a  computer  proj'ram 
written  for  the  IBM  709^  to  generate  raultigroup  cross- 
sections  using  the  or  B^  approximation.  It  will  cal¬ 
culate  these  constants  for  as  many  as  32  fast  groups. 
Available  as  output  are  diffusion  coefficients,  a  full 
scattering  matrix  which  includes  inelastic  and  (n-2n) 
scatter,  age,  and  constants  for  one  fast  group. 

There  are  three  options  for  calculating  cross-sec¬ 
tions;  one  may  use  the  equations,  use  the  B^  equations, 
or  road  in  a  flux  spectrum.  In  the  P^  approximation 
(Ref  1 4:336)  the  scattering  frequency  finctlon  and  flux" 
are  expanded  in  Legendre  polynomials,  and  the  P^  equations 
are  then  integrated  over  space  to  yield  two  energy  depen¬ 
dent  equations.  These  two  equations  are  then  solved  iter¬ 
atively  for  the  flux  and  current,  and  the  fine  (sixty- 
eight)  group  cross-sections  are  weighted  appropriately. 

The  B^  approximation  (Ref  9)  uses  a  Fourier  trans¬ 
form  of  the  one-dimensional  Boltzmann  equation,  with  the 
transformed  flux  and  frequency  function  expanded  in 
Legendre  polynomials.  Given  the  wave  number  "k" ,  the 
transformed  equations  for  flux  and  current  are  solved 
iteratively.  This  method  is  the  most  accurate  for  homo¬ 
genous  materials.  It  assumes  that  the  energy  spectxnim  of 
the  source  term  is  Independent  of  position.  To  get  some 
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Idea  of  the  dependence  of  the  approximation  on  the  wave 
number  used,  k  was  taken  as  the  square-root  of  the  trans- 

2  *4* 

verse  buckling  (B  and  then  of  the  total  buckling 

(B  )^.  The  change  in  cross-sections  for  water  was  very 

V 


slight . 

The  16-group  cross-sections  given  by  Yiftah  et  aJ^ 
(Ref  17)  wore  used  for  comparison  with  the  GAM-I  values; 
both  use  the  Cranberg  fission  spectrum  for  U-235  (Ref  ?)• 
Since  GAM-I  will  accept  any  fast  energy  groups  (up  to  32), 
Yiftah* s  group  structure  was  used  by  the  program  to  gen¬ 
erate  a  set  of  cross-sections.  Since  Yiftah  uses  a  dif¬ 
ferent  flux  spectrum,  the  results  were  not  expected  to 
be  identical.  The  two  sets  were  in  fairly  good  agreement, 
however,  especially  at  very  high  energies;  hence,  more 
faith  could  be  put  in  the  GAM-I  values,  and  this  code  was 
used  exclusively. 
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Appendix  C 

Personal  Correspondence 

The  two  letters  in  this  appendix  were  received  from 
workers  at  Los  Alamos  Scientific  Laboratory,  New  Mexico. 
They  are  self-explanatory  in  their  content,  the  first  re- 
ferrin^r  to  Ref  13  and  the  second  to  Ref  2,  as  applied  to 
the  derivation  of  S4C40;  both  are  valuable  as  sources  of 
additional  detail  for  the  above  references.  • 


61 


GNE/PHYS  65-10 


IN  M»LV 

■cm  to; 


UNIVERSITY  OF  CALIFORNIA 


LOS  ALAMOS  SCIENTIFIC  LABORATORY 
(CONTIUCT  W-7405-ENG-36) 

P.O.  Bo*  1663 

LOS  ALAMOS.  NEW  MEXICO  87341 


T-4 


January  4,  1965 


Lt.  Lee  H.  Livingston 
AFTT  Box  5368 
Wright -Patterson  AFB 
Ohio  45433 

Dear  Lt.  Livingston: 

In  reply  to  your  letter  of  December  2,  1964  I  would  like  to  offer  the 
following  comments: 

A  A 

Question  1:  Concerning  the  objects  f  and  3  of  Chapter  4.  On  page  75 

A  '  A 

(LA-2595)>  F  and  3  are  defined  by  Bqn.  4.8,  4.9  and  the  text  immediately  below. 
The  crucial  point  here  is  the  recognition  that  the  transport  equations  can  be 
transformed  into  an  "exact  differential"  form,  which,  when  integrated  over  the 
independent  variable,  depends  only  upon  the  end  points  (Condition  A).  How  is 
this  done?  If  one  uses  Bqn.  (4.7)  as  a  guide  in  writing  the  transport  equation 
terms  (see  Table  4.1),  then  the  application  of  Condition  A  (Bqn.  4.8)  becomes 
a  trivial  identity  when  the  integration  is  performed  over  the  applicable  inde¬ 
pendent  variable.  Bqn  (4.9)  is  likewise  an  identity.  However,  its  difference 
approximation  is  only  to  first  order. 

We  could  write,  equally  well,  Bqn  4.10,  as 

T(x)n(x)  -s  ^f(x)N(x)  =  f(x)  11+  N(x)  ~ 


a*  -j* 

ave 


-  V 


N 


+  N, 


3 
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where  f(>.)  is  some  analytic  function  occurring  in  the  transport  equation  and 

A  fflVP  *  /^f\ 

f  ^  8  3{^j  .  However,  we  have  still  to  determine  an  appropriate  form 

^Ax  ’  ^  Wave 

for  f  aiid  p  consistent  with  the  rest  of  the  approximation. 

The  piiilosophy  invoked  here  is  that  either  way  of  writing 

T(x)n(x)  =  ~  f(x)N(x)  =  f(x)  11+  N(x)  II 


in  difference  approximation  must  he  equivalent  in  its  application  to  the  re¬ 
mitting  eq\iations. 

In  Eqn.  (4.11),  the  integral  is  approximated  hy  the  averaged  value  of  the 
integrand  (us..ig  the  right  hand  side  of  Eqn.  (4.9))  with  an  unknown  (but  deter¬ 
minable)  weight  function  w.  Eqn.  (4.12)  performs  the  same  order  of  approximation 
but  with  the  alternate  form  (Condition  A).  Thus,  one  can  obtain  the  recursion 
relationship  given  by  Eqn.  (4.15). 

This  technique  is  just  that  required  to  guarantee  exact  conservation  of 
neutrons,  a  crucial  requirement  for  a  transport  code. 

Question  2;  Equivalence  of  Eq.  (4.32)  of  LA-2595  and  Eqn.  (?)  of  Carlson’s 
Numerical  Solutions  of  Neutron  Transport  Problems.  They  are  exactly  the  same 
equations  except  for  a  notational  change  in  the  angular  coordinates,  labeling 
of  the  interval  subscripts  and  use  of  the  diamond  differencing  assumption 
Eqn.  (4.l8).  We  now  indicate  the  details  for  clarification;  Write  Eqn.  (4.32) 
(ignoring  4th  and  5th  terms  for  the  moment)  as 


u 

V. 

1 


(A. -  A.N.)  + 
1+1  1+1  1  1 


1 

V. 


V.  ^\+l\+l 
1 


Eqn.  (7)  is  given  by 

(ti  cos  ^  +  “HA.  sin  ^,/a.  +  h.  +  u  A./\  )  N. 

'  ra  1  j  1  I  K  1 

+  (-ri  cos  cp^  +  sin  +  h^  +  jl  A^/a^^) 


-  2(TAj,  sin  (A^)  Nj_3_ 


_A. 


(2) 
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Eqn.  (1)  and  (2)  first  of  all  differ  in  the  angular  variable  notation,  so  let 
me  identify 


U  =  r|  cos  cp 


m 


ri  =  Ti  sin  CD, 


a'  =  h 


i=.jrT-2 


(3) 


Then  Eqn  (2)  becomes 


A.  A. 

2Ti'  N.  -  -  2§'  ~  N.  =  A.  S. 
A  j-l  A^  k-1  11 


(^) 


Next  divide  out  A.,  and  introduce  V.  =  v.  «  A.AA,  ,  yielding 

i  i  ijK  i  j  ic 

;4-A  A  +  S-AA  +  A-AA  +  ri  N.  +  {  -  JJ-A  A  +  ,4"^  A  +  N.  ^ 

i^V.  OK  V.  V.  1  j  A^/  1  V.  jk  V.  V..  1  j  A,/  i-1 


2 


-  ^  A  A  N.  ,  -  2  ~AA.  N.  ,  =  3. 

V,  i  k  j~l  1  0  k-1  1 


Now  regroup  the  terms  to  obtain 


^ ' 


‘A  ■  ”i-i>  *  <”1  ^  “i-i  - 


^  * 


+  (N.  +  N.  ,  -  2N,  J  +  C. 

V.  X  j  1  ,i-l  k-1  i 


C-.  (’iiy  =  s. 


(6) 


(5) 
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If  we  utilize  the  diamond  difference  condition  Eq.n.  (4.25) 


2N.'  =  N.  +  =  N.  +  =  N.  +  .  • 

1  1  1+1  j  J+1  t  L+1  m  JW-1 


in  the  above  equation,  then 


N.  +  N.  =  N,  +  N,  ,  =>  N.  +  N.  ,  -  2N,  =  N.  -  N, 

1  1-1  j  j-1  1  1-1  j-1  j  j- 


j  J-1 


«i "  “i-i  =  "k "  ”k-i  "i "  "i-i  -  ®k-i  “  \  -  “k-i 


N  +  N 

i 

2  ^ 


and  (6)  becomes 


Vj(“k  -  \.i>  *  Vi  -  Si  (9) 

Now  identify  the  area  factors  ^  1  ~  ^i^^  1  ~ 

A. A.  and  (9)  becomes 

^  o 

r  -  \-i  "i-i’ "  r:  <s/j  -  Sj-i  "i-i’  *  r '  Vk  -  “k-i  “k-i* "  "i"!  =  Si 


which  is  just  Eqn.  (l)  offset  by  one  interval.  Notice  that  Carlson  iises 

A.  =  r.  -  r.  ,  whereas  LA-2595  uses  A.  =  r_,  -  r,. 

1  1  1-1  1  1+1  1 


Question  5:  IVhy  do  the  fourth  and  fifth  terms  of  Eqn.  (4.52)  drop  out? 
By  definition  for  x,  y,  2  geometry,  Eqn.  (4.55) ^ 
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=  AytiiZ  ^  function  of  x 
^  function  of  y 
C,  =  ^  function  of  z 


The  subscript  convention  announced  was  i€x,  jey,  kez.  Thus,  A^  =  for 

all  i,  B  =  B...  for  all  C.  =  C.  for  all  k.  The  cross  section  perpen- 
j  j+l  k  K+l 

dicular  to  a  given  coordinate  is  independent  of  that  coordinate  since  the 
local  geometry  is  flat.  This  is  true  only  for  (x,y,z),  since  no  curvature 
terms  are  present  in  the  operators. 

The  condition  on 


Question  4;  Concerning  values  of  to  be  used. 


W,  is 
■Cm 


•tm 


=  1  when  siaaned  over  the  lowest  order  symmetry  of  the  sphere 


used  (octant,  hemisphere  etc.).  In  Table  4.5  this  is  two  octants,  and  hence 


^  W,  is  used.  If  a  hemisphere  is  necessary 


If  the  full 

sphere  (x,y,z)  is  used,  j  is  the  appropriate  quantity.  The  Appendix  A 
may  be  of  some  help.  Be  careful  of  direction  cosines  when  taking  the  full 


spnere . 

Finally  some  general  remarks.  Your  proposed  project  is  very  ambitious. 
You  should  realize  that  about  a  factor  25  increase  in  running  time  is  involved 
in  going  from  rz  to  xyz  geometry  due  to  increased  number  of  directions  and 
additional  required  space  points. 

We  include  an  rz  Fortran  Code  listing  for  your  information  and  the 
associated  flow  diagrams.  Tapes  and  cards  can  be  made  available  if  that 
would  help. 

If  we  may  be  of  further  assistance,  please  let  us  know. 


CEL/ml 

Enc;  Fortran  code  listing 
Flow  diagrams 


Sincerely  yours, 


Clarence  E.  Lee 
T-Division 


1 
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0?  REPLY 
REFER  TO: 


T-1 


March  1,  I965 


Mr.  Lee  H.  Livingston 
5466  Access  Road 
Dayton,  Ohio  -  4543I 

Dear  Itp.  Livingston: 

Bengt  Carlson  has  asked  me  to  answer  ycur  questions  of  I8  February. 
FrcHn  your  letter  I  gather  that  you  are'  referring  to  the  paper  in  Volume 
XI  of  Proceedings  of  Syngposia  in  Applied  Mathematics.  The  selection  of 
the  discrete  angles  (p  is  described  in  paragraph  5  onjpage  224  of  the 
proceedings.  In  equation  6  on  page  228  the  angles  qi^  and  9^  are  both 
the  same  angles,  i.e.,  The  absence  of  the  subscript  is  a  conven¬ 

tion.  The  orientation  of  these  angles  is  sketched  for  n=4  on  the  en¬ 
closed  page. 

■The  niaterial  described  in  the  above  reference  is  beidly  out  of  date. 
For  a  description  of  methods  in  use  in  1964  I  suggest  you  study  Los 
Alamos-  Scientific  Laboratory  Report  LA-2996.  Angular  orientations  and 
direction  sets  are  described  in  LA-3186,  a  copy  of  which  is  enclosed. 

I 

If  you  have  more  questions,  please  feel  free  to  write  to  me. 


Sincerely  yours. 


Dr.  K.  D.  Lathrop 


KDL:sv 

end:  1  graph  +  IA-3186 
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